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Abstract 

The adhesion approximation is a simple analytical model suggested for explanation of the 
major geometrical features of the observed structure in the galaxy distribution on scales from 
1 to (a few) x 100 Mpc/h. It is based on Burgers' equation and therefore allows analysis in 
considerable detail. A particular version of the model that assumes the infinitesimal viscosity 
naturally results in irregular tessellation of the universe. Generic elements of the tessellation: 
vertices, edges, faces and three-dimensional tiles can be associated with astronomical objects 
of different kinds: clusters, superclusters and voids of galaxies. Point-like vertices contain the 
most of the mass and one-dimensional edges (filaments) are the second massive elements. The 
least massive are the two-dimensional faces and tiles (voids). The evolution of the large-scale 
structure can be viewed as a continuous process that transports mass predominantly from the 
high- to low-dimensional elements of the tessellation. For instance, the mass from the cells flows 
into faces, edges and vertices, in turn the mass from faces flows into edges and vertices, etc. At 
the same time, the elements of the tessellation themselves are in continuous motion resulting in 
mergers of some vertices, growth of some tiles and shrinking and disappearance of the others as 
well as other metamorphoses. 

1 Introduction 

We live in an expanding universe theoretically predicted by Russian mathematician A. Friedman in 
1922 [14] and independently discovered by American astronomer E. Hubble in 1929 [15]. It means 
that in the perfectly homogeneous universe the relative velocity, v , between any two particles or 
objects would be described by the Hubble law 

v = H(t)r, (1) 

where r is the distance between the objects, the velocity is directed along the straight line connecting 
the objects. In the real inhomogeneous universe this law can be applied only at sufficiently large 
distances. A positive function Hit) is called the Hubble parameter, it characterizes the rate of 
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expansion of the universe. Its present value H = H(t ) 71 km/s/Mpc (Mpc stands for megaparsec, 
1 Mpc = 10 6 pc m 3,260,000 light years). Until the end of the millennium most cosmologists 
believed that the expansion of the universe was faster in the past and has been monotonically 
slowing down with time (i.e., H < at t > 0) due to the pulling effect of gravity. However the 
study of extremely distant supernovae stars revealed that the expansion of the universe recently (by 
cosmological standards) started to accelerate [31]. 

The acceleration of the universe can be easily explained if the most of mass in the universe is in the 
form of dark energy - a hypothetical form of energy uniformly filling all the space and experiencing 
a repulsive gravitational force. Although the physical nature of dark energy remains a mystery it 
is currently the most popular way of explaining the acceleration of the expansion of the universe. 
According to current measurements of the cosmological parameters dark energy accounts for 74% of 
the total mass-energy, 22% are due to dark matter - another puzzling component, while the ordinary 
matter often dubbed as baryonic matter makes only about 4% of the total mass. Dark matter is 
believed to consist of weakly interacting massive particles that interact with ordinary matter and 
themselves via gravity only for most intents and purposes. Although the exact kind of particles has 
been neither identified nor detected in physical experiments the hypothesis of dark matter is well 
founded. The gravitational interaction of dark matter is similar to that of ordinary matter. Dark 
matter plays a crucial role in the formation of the structure in the universe due to its dominance in 
density and therefore in generating attractive gravitational force. The universe is also permeated by 
the Cosmic Microwave Background (CMB) radiation that brings unique information about tiny initial 
perturbations of the uniformity of the universe that served as the seeds for the presently observed 
structure in the universe. However CMB radiation makes less than 0.01% of the total mass of the 
universe at present and has played a very limited role in the structure formation after the universe 
was roughly one thousandth of its present size. It has been also established that the geometry of the 
universe is very close to flat. The present universe is highly inhomogeneous on scales up to several 
hundred of Mpc as shown in Fig. 1. However, this figure does not show a true distribution of galaxies 
in physical space because the line of sight distances have been derived from the receding velocities 
via the Hubble law (Eq. 1). In the inhomogeneous universe the actual relation between the line of 
sight distance, r, and recession velocity, v b s , is given by the equation 

Vobs = H r + v p , (2) 

that involves an additional term, v p , called the peculiar velocity. The peculiar velocities are due to 
growth of density inhomogeneities in the universe that obviously requires a transport of mass from 
one place to another and therefore demands additional (with respect to uniform expansion) velocities. 
The recession velocities of galaxies, v Q b s , can be measured by making use of the Doppler effect, while 
the peculiar velocities cannot be easily measured. Therefore, the line of sight distance derived from 
the observed velocity r est = v o b s /H = r + v v j H gives only a crude estimate of the true distances of 
galaxies. 

The density distribution in physical space is better illustrated (albeit only in statistical sense) 
by the N-body simulations of the structure formation in the currently popular cosmological model 
referred to as the Lambda Cold Dark Matter (ACDM) model shown in Fig. 2. In cosmology the 
term large-scale structure is commonly referred to a distribution of galaxies on the scales roughly 
from 1 to a few hundred Mpc. Galaxies can not probe smaller scales because of discreteness, and 
on larger scales the galaxy distribution is monotonically approaching homogeneity with the growth 
of the scale. The redshift surveys reveal spectacular abundance of structures often described as 
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Figure 1: Sloan Digital Sky Survey Great Wall (top) compared to CfA2 Great Wall (bottom) at the same 
scale. The scale is shown in the velocity units (km/s), the distance in Mpc can be obtained from the Hubble 
law (Eq. 1), e.g., 10,000 km/s correspond to 141 Mpc. Each point represents a galaxy Adopted from [11]. 

filamentary, network, or bubble structure [10, 8]. The origin of the large-scale structure is one of the 
most important problems in modern cosmology. Many fundamental issues in physics, cosmology and 
astronomy ranging from speculations on the physical nature of dark matter, to the measurement of 
angular anisotropy of the microwave background radiation as well as establishing the epoch of galaxy 
formation join together here, see e.g., [30, 34]. 

Modern theory explains the formation of the structure in the universe as a result of growth 
of primeval tiny density fluctuations in the process known as gravitational instability. Primeval 
perturbations are assumed to arise as vacuum fluctuations during the very early stage when the 
universe was expanding exponentially (inflationary universe), see e.g., [28]. Afterwards, the density 
perturbations had a long evolution before they become galaxies, clusters of galaxies, superclusters 
and voids of galaxies. The formation of galaxies is a very difficult problem on its own. Many complex 
physical processes like star formation and supernova explosions are very important for understanding 
the galaxy formation. We shall discuss the mass distribution only assuming that galaxies are fairly 
good (though not perfect) tracers of mass on large scales. 

As long as the density perturbations are small by amplitude their evolution is described by the 
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Figure 2: The dark matter density field in a slab of thickness 15 Mpc/h obtained in the Millennium 
simulation [16] of the structure formation in the universe. 

linear theory of gravitational instability (see e.g. [30, 34]). The linear theory is quite simple and 
therefore is well understood. In particular, it predicts the rate of growth of the perturbations, D(t): 
5 = Sp/p oc D(t), where 5p = p — p and < 5 2 >^ 2 <C 1. However, in commonly used Eulerian form 
it is not very useful at the nonlinear stage when the amplitude of the density perturbations becomes 
large and the observable structures (sheets, filaments, and clusters of galaxies) form. 

The analysis of the evolution of perturbations at the nonlinear stage < S 2 > 1//2 > 1 becomes 
quite difficult. The most straightforward way to address the complexity of nonlinear evolution is 
to carry out three-dimensional N-body simulations (Fig.2). Usually in the simulations of this type 
the medium is assumed to consist of collisionless particles, in agreement with the hypothesis of dark 
matter in the form of weakly interacting particles. The trajectory of each particle in the simulation 
is numerically integrated in the gravitational field generated by all the rest particles. Boundary 
conditions are commonly assumed to be periodic that imitate the infinite universe. 

Here we outline another approach to the problem of the large-scale structure in the universe. 
We present two approximate analytic solutions to the set of PDEs describing the growth of density 
inhomogeneities in an expanding universe. One suggested by Zel'dovich in 1970 [38] is known as the 
Zel'dovich approximation and the other is based on Burgers' equation [5] and is often referred to 
as the adhesion model [12, 13] (see also [34]). A particular version of the adhesion model naturally 
describes the structure in the universe as an irregular tessellation or quazi-Voronoi tessellation [27]. 
The generic elements of the tessellation in 3D (vertices, edges, faces, and tiles) can be associated 
with the observed structures in the three-dimensional distribution of galaxies (clusters, superclusters 
and voids of galaxies). 
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Both the N-body simulations and approximations require initial conditions. In the linear regime 
the density fluctuations are assumed to be a realization of a Gaussian random field specified by the 
power spectrum and the amplitude. The current measurements of the angular fluctuations in the 
temperature of the microwave background radiation by WMAP (Wilkinson Microwave Anisotropy 
Probe) are in excellent agreement with the assumption of Gaussianity. [22]. The amplitude of the 
temperature fluctuations suggests that the scale where the density fluctuations have recently reached 
nonlinearity is about 3-6 Mpc depending on exact definition of the scale of nonlinearity. This 
is also in a good agreement with the observations of the large-scale distribution of galaxies. The 
shape of the power spectrum of the initial (i.e., linear) density fluctuations is fully determined by 
the parameters of the cosmological model. At large scales (small k) it behaves as P(k) oc k while at 
very small scales it is approximately P(k) oc k~ 3 with smooth transition between these limits. It is 
approximately P(k) oc A; -1 ' 5 at the present scale of nonlinearity, see the discussion bellow and Fig.5. 

Here we discuss the models that deal only with the dominating component of the mass able to 
cluster i.e., the dark matter. The luminous objects in the universe are made of baryons, therefore 
it is important to understand the dynamics of the baryonic component of mass as well. However, 
this problem is much more complicated and is beyond the scope of this review. We refer the inter- 
ested reader to paper [17] that offers an interesting generalization of the models described here and 
addresses the dynamics of baryons. 

The rest of the paper is organized as follows. The equation describing the evolution of density 
inhomogeneities are present in Sec. 2. Sections. 3 and 4 briefly describe the Zel'dovich approximation 
and the adhesion model respectively. Section 5 provides a short summary. 



2 Basic Equations 

The evolution of density inhomogeneities can be described by a system of three partial differential 
equations comprising the continuity, Euler, and the Poisson equation (see e.g. [30, 34]). In order to 
exclude the uniform expansion of the universe comoving coordinates, x, and peculiar velocity, u p are 
commonly used. They are defined by the following relations 

r = a(t)x, r = H(t)r + u p . (3) 

A monotonically growing function a(t) called the scale factor describes the uniform expansion of 
the universe, and H(t) = dlna/dt is the Hubble parameter introduced by Eq.l. The scale factor 
is completely determined by the parameters of the cosmological model. For instance, a(t) oc t 2 ^ 3 in 
the Einstein - de Sitter universe where A = 0, and Q = 1, while in more realistic ACDM model it is 
a complicated but known function of time [6]. In the course of evolution neither peculiar velocities 
no gravitational potential reach the relativistic values, therefore the use of classical mechanics and 
Newton's gravity law is perfectly appropriate. In terms of the comoving coordinates and peculiar 
velocities the equations describing the gravitational instability in the expanding universe are as 
follows: 

the continuity equation 

+ I V • (p Up ) = -3Hp, (4) 

the Euler equation 

+ -(up • V)u p = — V0 - Hu p , (5) 
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and the Poisson equation 



lv 2 = 47rG(p-p), (6) 



where p and p are respectively the density and mean mass density; (J) is the gravitational potential 
due to the inhomogeneities of density; G is the gravitational constant. 

The pressure is neglected since we study the medium interacting only gravitationally. Additional 
terms on the right hand side of the continuity and Euler equations (i.e., —3Hp and —Hu p respectively) 
are due to the expansion of the universe and the factor 1/a is due to differentiation with respect to 
the comoving coordinates x: V = d/dxi = a ■ d/dri. The initial conditions are small density and 
smooth velocity perturbations imposed on a homogeneous density distribution. 

As long as the amplitude of the density perturbations remains small their evolution can be 
analyzed in the linear approximation obtained by the linearization of the above equations. The exact 
solution of the linearized system has one growing mode which is the major object of our analysis 
and two decaying modes that can be ignored. The velocity in the growing mode is a potential vector 
field which is proportional to the gradient of the linear gravitational potential: vi in oc — V0u n . In 
the linear regime the spatial structure of the perturbations (in the comoving coordinates) remains 
unchanged and its amplitude is proportional to the growing solution 5 oc D{t) also determined by 
the cosmological model [6]. For example, D{t) oc a(t) oc t 2 ^ 3 in the Einstein - de Sitter model, in 
general D(t) is a monotonically growing function of time. 

Equations 4-6 describing the evolution of inhomogeneities in the universe become more convenient 
for further analysis after the following transformation of variables are made 



p = a 3 r], u p = aDv, <p = (^Q d 2 D^j if. (7) 
In addition, we shall use D instead of time, t. The resulting set of equations becomes 

dD dx; 



dvi dv 

dD +Vk dx k 2Df 




^+ v k^T = -oTT^I— ( 8 ) 



d 2 ip 6_ 
dxf ~ D' 

Here Qo = p m j p c at present time (it is the ratio of the total mean mass density to the critical value 
p c = 3Hq/8ttG), f = din D/ din a and S — (rj — fj)/fj where fj = flop c = const; we also assume the 
summation over repeated indices. 

The second equation in set (8) becomes even simpler after introducing the total derivatives 
d/dD = d/dD + v k d/dx k 

dvi 3 n ( dip \ 

dD = -2Dp\^ + V T () 
One can easily check that in the linear regime (i.e., while < 5 2 > 1 ^ 2 ^C 1) the only growing solution 
to the set of equations (8) is given by 

K*H - 

v(q,£>) = v (q) = -V,$ (q), (10) 
(f(q,D) = $ (q), 
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where q is the Lagrangian coordinates of a fluid element. The solution (10) gives the density contrast, 
velocity and potential in terms of the Lagrangian coordinates. In order to find these fields in the 
Eulerian space one needs to solve the equation of trajectories 



x(q,£>) = q + J D-v (q) 



(11) 



for q = q(x, D) and substitute it in solution (10). The field $o(q) is determined by the perturbation 
of density in the linear regime which is assumed to be a realization of a Gaussian random field. The 
power spectrum of the field that solely defines all its statistical properties is derived from observations 
of the fluctuations of temperature in the Cosmic Microwave Background radiation (see e.g., [28]); it 
is shown in Fig. 5. 

Although velocity Vi remains constant for each particle in the linear regime the physical velocity 
Up varies with time (Eq.7). In the growing solution (Eq.10) gravitational force proportional to 
—d(p/dxi, exactly balances the drag force due to the uniform expansion of the universe Hu p oc v 
and thus the right hand side of Eq.9 vanishes to linear order. 

3 The Zel'dovich Approximation 

Zel'dovich (1970) derived the Lagrangian form of the linear theory and extrapolated it to the be- 
ginning of the nonlinear regime [38] (see also [34]). This extrapolation assumes that the initial 
perturbations are smooth in the sense that the initial power spectrum must decrease faster than k~ 3 
at k — > 0. In realistic cosmological models e.g., in the ACDM model the initial power spectrum 
P(k) oc \n(k)k~ 3 up to very large k corresponding to cosmologically irrelevant scales and thus the 
Zel'dovich approximation must be properly modified in order to be useful on the scale of the clusters 
and superclusers of galaxies. This will be discussed later. 

The Zel'dovich approximation can be conveniently formulated as a mapping from Lagrangian 
space L{q} to Eulerian space E{x) defined by Eq.ll. A two-dimensional illustration of the mapping 
is shown at four stages of the evolution in Fig.3. The stages are marked by a =< 5f in >^ 2 where 
5 = [p(x, t) — p(t)]/p(t) = [?y(x, D(t)) — fj]/fj is the density contrast. By convention, the parameter a 
is computed in the linear theory and therefore a oc D(t). 

Utilizing conservation of mass i]d 3 x = fjd 3 q Zel'dovich has derived the density as a function of 
D(t) and the Lagrangian coordinate q 



where Ai(q), A2(q) and Aa(q) are the eigenvalues of the deformation tensor dij = d 2 $o/dqidqj. In 
cosmology the initial condition is usually characterized by the spectrum P$(k) (Fig. 5); it is obviously 
related to the power spectrum of the potential as P$ = k 4 P^ . Similarly to the linear stage one can 
find the density distribution in the Eulerian space by solving Eq.ll for q and substituting it in Eq.12. 
For realistic initial conditions it requires numerical calculations. 

It follows from Eq.12 that the first objects are formed around the peaks of the largest eigenvalue 
(we assume that A; are ordered at every q: Ai > A2 and A2 > A3) and have very oblate shapes 
because generic peaks of Ai always have two other eigen values different from Ai. These objects are 
known in cosmology as "Zel'dovich's pancakes". This prediction of the Zel'dovich approximation 
differs from the extrapolation of the Eulerian linear theory that predicts the formation of the first 
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[1 - £>(f)Ai(q)][l - £>(f)A 2 (q)][l - D(t)X 3 (q)\ ' 
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Figure 3: The distributions of particles at four stages of evolution described by the Zel'dovich approximation 
in two dimensions. The panels are marked by a =< Sf in > 1 / 2 , the variance of the density contrast evaluated 
to the linear order, so a oc D(t). 

objects from the peaks of 5 = Ai + A 2 + A 3 . In practice the difference in most cases is not large 
due to obvious correlation between 5 and Ai that causes the peaks of 5 be spatially close to the 
peaks of Ai. Three-dimensional gravitational N-body simulations are in a good agreement with the 
prediction of the Zel'dovich approximation [35]. Pancakes in the collisionless dark matter originate 
as the three-stream flow regions bounded by caustics, the surfaces of formally infinite density. The 
shape and other characteristics of the pancakes as well as of other generic types of structures (Fig. 3) 
are determined by catastrophe theory [2]. 

The Zel'dovich approximation proved to be very good until orbit crossing when caustics form and 
the multi-stream flows occur i.e., up to the stage corresponding to a m 1 (see e.g. [34] and references 
therein). At a > 1 Eq.ll predicts broadening of the multi-stream flow regions considerably faster 
than it happens in N-body simulations [9, 18, 29] as can be seen in Fig. 4. Compare two panels on 
the right in Fig.3 with two panels of Fig. 4. Despite the difference in the realizations of the initial 
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Figure 4: The distribution of particles at two stages <r = lon the left and a = 2 on the right obtained in a 
two-dimensional N-body simulation. Adopted from [32]. 

conditions one can see that the correct gravitational gravitational force in the N-body simulations 
results in much thinner pancakes than predicted by the Zel'dovich approximation. This discrepancy 
motivated the development of the adhesion model that will be discussed in the following section. 
Here we briefly describe a simple modification that allows to utilize the Zel'dovich approximation in 
the case when the small scale power cannot be ignored [7, 25] . 

N-body simulations both in 2D and 3D have shown that the largest structures depend very little 
on small-scale perturbations [3, 23, 24]. The most important factor that affect the shapes and general 
appearance of the large-scale structure are the initial waves just reached nonlinear regime, i.e., on 
scales where rms fluctuations are just becoming nonlinear [23]. This observation allows to use an 
auxiliary model that yields the large-scale structures similar to the original model but has no small- 
scale structures [7, 25]. The scale that separates large-scale structures from small-scale structure is 
the scale of nonlinearity. 

For instance, the initial power spectrum in the currently popular ACDM cosmological model is 
shown in Fig. 5 by the solid line while the spectrum filtered at the present scale of nonlinearity is 
shown by the dotted line. If the realization of amplitudes and phases are same in both models then 
the auxiliary model with filtered power spectrum would have a very similar large-scale structure to 
one in the model with unfiltered power spectrum. The dotted line in Fig. 5 shows the spectrum filtered 
with a Gaussian window; a simpler version of filtering when the power spectrum set to zero on scales 
with k > k n \ gives a similar result. The auxiliary model with filtered power spectrum provides an 
option of using the Zel'dovich approximation and therefore better understanding of complex nonlinear 
processes that control the structure formation. The price of this simplification consists in missing 
all structures on smaller scales that are abundant in the model with unfiltered initial spectrum. The 
adhesion model is an attempt to remedy this drawback of the Zel'dovich approximation. 

Finally, it is worth noting that both the Zel'dovich approximation (Fig. 3, top right panel) and 
N-body simulations (Fig. 4) suggest that the universe can be viewed as an irregular tessellation where 
'tiles' are represented by voids, the regions of low mass and/or galaxy densities, where the flow 



9 




LCDM 
Filtered LCDM 



0.0001 



1e-05 0.0001 0.001 0.01 



0.1 



10 



100 



k (h/Mpc) 



Figure 5: The initial power spectrum in the ACDM cosmology is shown by the solid line, the spectrum 
filtered at the current scale of nonlinearity is shown by the dotted line. 

consists of only one stream. 



The general idea of the adhesion model is very simple. We wish to use the Zel'dovich approximation 
everywhere except the regions of multi-stream flows. By adding a diffusion term to the Euler equation 
one can avoid the orbit crossing and therefore suppress the formation of the multi-stream flow regions. 
The advantage of this modification is the suppression of too fast growth of the thickness of pancakes 
observed in the Zel'dovich approximation but it is achieved at the coast of noticeable change of the 
dynamics inside multistream flow regions. A small coefficient of viscosity guarantees agreement with 
the Zel'dovich approximation before the orbit crossing. 

Since the motion is potential one can introduce the velocity potential v = — V$. Then assuming 
that the gravitational potential is approximately equal to the velocity potential ip ~ $ and also 
adding a viscosity term zA7 2 v to the right hand side of the second equation in set 8 one obtains the 
equation of nonlinear diffusion [12, 13] 



Generally speaking the viscosity term needs not to be in the form of Eq.13 but choosing this particular 
form one obtains Burgers' equation that has an exact analytic solution [5] . For the potential motion 
v = — V$ Eq.13 can be solved by performing the Hopf-Cole substitution $(x, D) = —2v log U (x, D). 
As a result Eq.13 adopts the form of the linear diffusion equation 



4 The Adhesion Model 



<9v 
dD 



+ ( v • V)v = zA7 2 v. 



(13) 




(14) 
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Solving Eq.14 for the velocity one obtains 

j^ 3 g [(x-q)/D] exp[S(x,D;q)/2»/] 
nX '^ J exp[5(x,L>;q)/2H ' 

where the "action" 

S(x,D;q) = $ ( q )-^. (16) 




m < 15.5 
10258 galaxies 

Figure 6: Simulation of galaxy distribution in redshift space using the adhesion approximation with small 
but finite viscosity. Adopted from [37]. 



In cosmology the adhesion model has been used in two forms: one assumes a small but finite 
value of the viscosity parameter v and the other assumes it is infinitesimal. For small but finite v the 
trajectory of a particle can be determined by solving the integral equation and the resulting density 
can be determined from the continuity equation [36, 37, 29, 26]. An example of the structure obtained 
by this method is shown in Fig.6 where the authors [37] plotted the simulated galaxy distribution in 
a slice through the simulation box similar to the plots obtained from galaxy redshift catalogs as in 
Fig.l. An interesting attempt to explain the magnitude of the coefficient of viscosity from dynamics 
in multistream flow regions was made in [4]. 

We will discuss a different approach that assumes the infinitesimal value of the viscosity parameter 
v — > in more detail because it has direct connection to tessellation and tiling. In this case the 
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Figure 7: The illustration of the geometrical technique that finds the Eulerian position of a particle using 
the initial potential <J>o(q)- For a given value of D the paraboloid (Eq.17) descends on the surface of the 
initial potential <3?o until it is tangent to the surface at point q by adjusting constant Pq. The projection of 
the apex of the paraboloid denotes the Eulerian position of the particle. Adopted from [20]. 

integrals in Eq.15 can be evaluated using the method of steepest descents [12, 13, 20, 21]. The result 
has a remarkable geometrical interpretation illustrated by Fig. 7 in two-dimensional case. The initial 
Lagrangian coordinates q p and the current Eulerian coordinates x p of a particle at a chosen stage 
denoted by D, can be related simply by descending a paraboloid 



with apex at x = x p and adjusting correspondently the constant Pq until it is in contact with the 
surface of the initial potential $o(q) & t some point. The coordinates of the contact point are the 
Lagrangian coordinates of the particle q p . This procedure must satisfy an additional constraint: 
the paraboloid must not cross $ a t any point. Thus, the rule is: the tangential contacts of the 
paraboloid with the potential are welcome but crossings are strictly forbidden. At the linear stage 
D is small and the curvature of the paraboloid is greater than the typical curvature of $o, thus this 
condition can be easily fulfilled. The Zel'dovich approximation is universally valid at this stage and 
the velocity filed is given by the following equation 



P(x,D;q) 




(17) 



D 



o- 



v(x,L>) 




(18) 
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Figure 8: An illustration of the geometrical prescription of descending a parabola onto the initial velocity 
potential in order to find the Eulerian positions of particles and knots in one-dimenaional case. The particle 
having Lagrangian coordinate qo in the uppermost panel has the Eulerian coordinate of the paraboloid apex. 
In the middle panel corresponding to a later stage the whole region between two points where the parabola 
is tangent to the initial potential is stuck into the knot m\ moving with the velocity shown by the arrow, 
its mass is m\ = fj(q2 — Qi)- The knot 777-2 is determined by the second paraboloid in the middle panel. The 
lower panel shows the knot M formed as a result of merging mi, 7772 and the mass between them. Adopted 
from [32]. 

where q(x, D) is the coordinate of the absolute maximum of the action S(x, D; q) at given x and D. 

As time passes the curvature of the paraboloid decreases, so it becomes wider. As a result 
the points q that are not accessible without violation of the crossing rule as can be seen in one- 
dimensional illustration in Fig.8. Both parabolas shown in the middle panel cannot descend further 
without crossing the initial potential <3> . What has happened to the corresponding particles? The 
particles between two contact points i.e., between qi and q 2 for one parabola and between q 3 and 
g 4 for the other have been stuck into point masses mi and m 2 respectively due to orbit crossing. 
Their coordinates are the coordinates of the apices of the corresponding parabolas and the masses 
are m 1 = 77(^2 _ Qi) and m 2 = 77(^4 — ^3). At later times the discrete masses formed earlier can grow 
by merging as well as by accreting the continuous medium between them. This is illustrated by the 
bottom panel. It is not hard to see that the transition from the state shown in the middle panel 
to one shown in the bottom panel requires the state when one parabola has three contacts with the 
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Figure 9: The left hand side panel shows the distribution of particles obtained in two-dimensional N-body 
simulation. The panel on the right shows the tessellation computed in the adhesion model from the same 
initial conditions. The tessellation is superimposed on the particle distribution. Adopted from [21]. 

potential simultaneously. At the next instant of time it becomes a little wider and looses the contact 
with the central peak. 

In two dimensions the paraboloid can be tangent to the initial potential $o simultaneously in 
one, two or three points. The first type corresponds to the particles that form continuous medium. 
The points of the second type have been stuck into lines shown by thick points in Fig. 7. These lines, 
the edges of tiles, are formed by the set of the apices of the paraboloids that are in contact with the 
initial potential at two points simultaneously. Finally, when the paraboloid is in contact with the 
surface of $o & t three point simultaneously its apex marks the knot (vertex) where the lines (edges) 
meet. The edges and vertices form the boundaries of irregular tiles of tessellation in two dimensions. 
In the course of the evolution some tiles expand while the other collapse and cease to exist. This 
goes hand in hand with merging of the knots. At the instant of a merger the paraboloid can be in 
contact with the initial potential at more than three points simultaneously [1]. 

Figure 9 demonstrates the agreement between the tessellation computed in the adhesion model 
and the N-body simulations in two-dimensional space. Figure 10 illustrates the effect of small-scale 
power. The initial potentials in both panels were generated with the same amplitudes and phases 
but the amplitudes for the short waves k > k cr were set to zeros in the case shown in the left hand 
side panel. Therefore, the initial potential in the left hand side panel is a smooth version of the 
potential in the right hand side panel. The both models evolved to the same stage. One can see that 
while the overall appearance of the tessellation is similar in both panels the there more small-scale 
details in the panel on the right. 

With a little imagination and help from the geometric model outlined above one can easily 
visualize the main features of the tessellation in three-dimensional space. The initial potential $o 
forms a three-dimensional hypersurface in four- dimensional space. A three-dimensional paraboloid 
(Eq.17) descends on the thre-dimensional hypersurface of the initial potential similar to one- and 
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Figure 10: Two panels show the tessellation of the plane by the adhesion model superimposed on the 
particle distribution obtained in 2D N-body simulation. The initial potential <I>o i n both panels have the 
same large-scale part but on the left hand side it does not have the small-scale part that is present in the 
potential on the right hand side. Adopted from [21]. 

two-dimensional cases. The rules are the same: it can touch the initial potential but must not cross 
it. Four generic situations are possible at a generic instant of time: it can have only one, two, three 
or four contact points simultaneously. If there is only one contact point then the corresponding 
particle is in free motion; the coordinates of the of contact point are the Lagrangian coordinates of 
the particle and the apex denores its Eulerian coordinates. In the remaining three cases the apex 
define a particle stuck on a surface, on an edge or at a vertex of the tessellation, if the number of 
contact points is two, three or four respectively. Again similarly to one- and two-dimensional cases 
at the critical instants of time when a tree-dimensional tile collapses and ceases to exist or two or 
more vertices merge the number of contact points can be greater than four [1]. 

The set of surfaces, filaments and knots define the tessellation of three-dimensioal space: volumi- 
nous low density tiles (voids) are separated by the faces (pancakes), the edges (filaments) are at the 
intersections of the faces, and the vertices (clusters of galaxies, halos in N-body simulations) are at 
the intersections of the edges. This tessellation can be viewed as the skeleton of the real structure. 

The both variants of the adhesion model (with small but finite and infinitesimal viscosity) have 
been tested against the gravitational N-body simulations in two and three dimensions. Both the 
N-body simulation and the adhesion model used the identical initial conditions and were compared 
at several stages of the evolution. The geometrical version of the adhesion model was tested against 
the two-dimensional N-body simulations with the initial power law spectra Ps(k) oc k n with spectral 
indices n = 2,0 and —2 and various cutoffs [21]. The version of the adhesion model utilizing a 
finite viscosity parameter v has been quantitatively compared to fully nonlinear, numerical three- 
dimensional gravitational N-body simulations in [26]. The results generally show a good agreement 
however they improve with the growth of the relative power on large scales. This is the case in the 
ACDM model as one can see from the initial power spectrum shown in Fig. 5. 
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5 Summary 
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Figure 11: Normalized initial gravitational potential at the centers of voids is plotted against the void 
diameter. Adopted from [32]. 

The adhesion model based on three-dimensional Burgers' equation (Eq.13) of the nonlinear dif- 
fusion has been found to provide a good explanation of the major large-scale features of the galaxy 
distribution in the universe. Its particular version assuming the infinitesimal coefficient of viscosity 
v — > provides a natural geometrical construction that can be characterized as " quasi- Voronoi tes- 
sellation" [27]. Comparison with gravitational two- and three-dimensional N-body simulations has 
shown a good agreement including advanced nonlinear stage. The adhesion model provides a natural 
qualitative explanation of the origin of the large-scale coherent structures such as superpancakes and 
superfilaments, as a result of coherent motion of clumps due to large-scale inhomogeneities in the 
initial gravitational potential [19]. The initial gravitational potential also affects the sizes of voids 
as shown in Fig. 11: the higher the peak of the initial potential the greater size of the voids formed 
around it. The formation and evolution of large-scale structure is described by the adhesion model as 
a two stage process [32]. During the first short stage matter falls predominantly onto pancakes and 
then moves along them towards filaments and then along filaments to the knots. At the end of the 
first stage the formation of the tessellation is completed and the bulk of the mass in the universe is 
located in knots (vertices) and filaments (edges) however a relatively small fraction of mass remains 
in pancakes (faces) and voids (three-dimensional tiles). The second stage is characterized by the 
deformation of the tiles and edges due to the dynamics caused by gravity . At this stage knots merge 
into more massive knots and small voids shrink and cease to exist giving space to growing large 
voids. Eventually almost all the mass concentrates in knots. Depending on the initial spectrum the 
knots themselves may move coherently in such a manner that they concentrate to superpancakes and 
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superfilaments. This situation is likely to occur in the ACDMmodel at the present time due to the 
shape of the initial power spectrum (Fig. 5). The superpancakes and superfilaments can be identified 
by applying the adhesion model to smoothed initial potential [21]. 

Acknowledgement The author is grateful to Rien van de Weygaert for many useful discussions 
related to the subject of this review. 
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